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Large-eddy simulation of conductive flows at low 
magnetic Reynolds number 

By B. Knaepen AND P. Moin 


1. Introduction 

In large-eddy simulations (LES), only the large-scale structures of the flow are simu- 
lated directly while the small-scale structures are taken into account through a model, 
referred to as the subgrid-scale (SGS) model. This results in significant reduction in com- 
putational requirements. The trade-off comes in the extra modelling effort that has to be 
produced in order to adequately take into account the discarded small scales structures. 
In the computational fluid mechanics community, the most widely used SGS model is 
the Smagorinsky eddy viscosity model (Smagorinsky 1963). This model has gained an 
even bigger practical interest following the work of Germano et al. (1991) in which the 
dynamic procedure was introduced. The dynamic procedure enables optimization “on the 
fly” of the arbitrary scaling factor that is inherently present in the original Smagorinsky 
model (more details below) and thus allows the model to automatically adapt to the flow 
being studied. 

In this paper we study the LES method with dynamic procedure in the context of 
conductive flows subject to an applied external magnetic field at low magnetic Reynolds 
number R rn . These kind of flows are encountered in many industrial applications. For 
example, in the steel industry, applied magnetic fields can be used to damp turbulence in 
the casting process. In nuclear fusion devices (Tokamaks), liquid-lithium flows are used as 
coolant blankets and interact with the surrounding magnetic field that drives and confines 
the fusion plasma. Also, in experimental facilities investigating the dynamo effect, the 
flow consists of liquid-sodium for which the Prandtl number and, as a consequence, the 
magnetic Reynolds number is low. 

Most of the previous works considering LES in the case of MHD flows have been 
directed towards flows at high R m , where the full non-linear MHD equations have to be 
used, or towards flows in the absence of an external magnetic field: Yoshizawa (1987); 
Zhou & Vahala (1991); Theobald et al. (1994); Agullo et al. (2001); Muller & Carati 
(2002). To our knowledge, the only attempt to study MHD turbulence at low magnetic 
Reynolds number from the LES point of view is due to Shimomura (1991). In that 
work, the author considers the case of magnetohydrodynamic turbulent channel flow and 
introduces a new SGS model designed to incorporate the effects of the applied uniform 
magnetic field. However, like the original Smagorinsky model, the model of Shimomura 
contains constants that have to be adjusted and are thus flow dependent. In the present 
paper, we show that it is possible to circumvent this problem by the use of the dynamic 
Smagorinsky model. 

Our attention is focused here on the case of homogeneous (initially isotropic) de- 
caying turbulence. The numerical simulations performed mimic the thought experiment 
described in Moffatt (1967) in which an initially homogeneous isotropic conductive flow 
is suddenly subjected to an applied magnetic field and freely decays without any forcing. 
Note that this flow was first studied numerically by Schumann (1976). It is well known 
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that in that case, extra damping of turbulence occurs due to the Joule effect and that 
the flow tends to become progressively independent of the coordinate along the direction 
of the magnetic field. Our comparison of filtered direct numerical simulation (DNS) pre- 
dictions and LES predictions show that the dynamic Smagorinsky model enables one to 
capture successfully the flow with LES, and that it automatically incorporates the effect 
of the magnetic field on the turbulence. 

Our paper is organized as follows. In the next section we summarize the LES approach 
in the case of MHD turbulence at low R m and recall the definition of the dynamic 
Smagorinsky model. In Sec. 3 we describe the parameters of the numerical experiments 
performed and the code used. Section 4 is devoted to the comparison of filtered DNS 
results and LES results. Conclusions are presented in Sec. 5. 


2. LES at low magnetic Reynolds number 

For homogeneous flows, the magnetic Reynolds number can be defined through the 
following relation: 

Rm = — , (2.1) 

n 

where u = ( UiUi ) /3 is the r.m.s. of the fluctuating velocity Uj, L is the integral length 

scale of the flow and is the magnetic diffusivity. R rn represents the relative importance 
of the non-linear terms and the diffusion term in the magnetic induction equation. In the 
limit of low R rn , the MHD equations can be simplified considerably Roberts (1967). It is 
indeed possible to close independently the momentum equation and take into account the 
effect of the magnetic field through an extra damping term. This approximation is known 
as the quasi-static (QS) approximation and along with the incompressibility condition, 
diUi = 0, reads, 

/ Tjext\2 

d t Ui = -di(p/p ) - UjdjUi ^ — A ~ 1 d z d z u i + vAu^ (2.2) 

where p is the sum of the kinematic and magnetic pressures, p is the fluid density, v the 
kinematic viscosity, B ext is the applied external magnetic field and A -1 is the inverse 
of the Laplacian operator. Note that B% xt has by convention been aligned with the z- 
direction and expressed in Alfven units. 

LES equations are obtained by filtering (2.2). The filtered velocity u, is defined by, 

Ui{x) = J G(x,y) uj(y)dy, (2.3) 

where G is a smoothing kernel that eliminates the small scale part of iq and satisfies the 
relation, 


J G(x, y)dy = 1. 

In terms of the filtered velocity, Eq. 2.2 can be written as, 

d t Ui = - 80/ p ) - UjdjUi 5 A ~ 1 d z d z u i + vA Ui - djfij , 


(2.4) 


(2.5) 



where. 
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Tij = UiUj - UiUj, (2.6) 

is the subgrid-scale stress tensor. In order to close (2.5), T i7 has to be expressed only in 
terms of the filtered velocity. 

Is is interesting to note that T i7 does not depend explicitly on the magnetic field. 
Indeed, the magnetic contribution, being a linear term in (2.2), does not require an 
explicit SGS counterpart in the LES equation. In the case of non-conductive flows, the 
most widely used model for Tij is the Smagorinsky model: 

= -2C s A 2 \S\Sij, Sjj = ^( dfxij + djuf), S = ^J~2S~S~j, (2.7) 


where A is the width of the filter G and C s is the Smagorinsky constant. However, as 
was observed by Shimonura in Shimomura (1991), model (2.7) with C s optimized for a 
non-conductive flow, is not adequate in the present context. This is easily understood 
when one refers to the fact that the applied external magnetic field has the tendency 
to suppress non-linear transfers in the velocity field. Thus, the effect of the SGS stress 
tensor needs to be decreased for conductive flows when the external magnetic field is 
switched on. This should result in a lower optimal value for the Smagorinsky constant. 

Contrary to what was done in Shimomura (1991), we will not explicitly modify the 
Smagorinsky model to incorporate the effect of the magnetic field but rather we will 
make use of the dynamic procedure (Germano et al. 1991) to optimize the value of C s . 

To define the dynamic procedure one introduces a second, coarser filter called the 
test-filter which we denote by ~. Applying this filter in addition to (2.5) yields, 


d t Ui = -di(p/p) - UjdjUi - 


-O- 2 Hr x d z d z Ui + z'AtZj - djTij - djLij, (2.8) 


where L,j = UiUj — UiUj is the Leonard tensor which does not require any modeling since 
it is expressed in closed form using the filtered velocity Wj. The sum Tij + Lij represents 

the SGS stress tensor Ty of the combined 7TT filter and we have the well-known Germano 
identity, 


t^. (2-9) 

Assuming that Tii and Hi are self-similar, a suitable model for T. t j should be, 

= -2G s A 2 |S|§y, + djTii), 13 = \fl (2.10) 

where A is the width of the 7TT filter. When one models r,; ? and T,j using (2.7) and 
(2.10), the Germano identity will undoubtedly be violated. However, the constant C s 
can be chosen in such a way as to minimize the difference (in the least square sense) 
between both sides of (2.9). Assuming homogeneity in all directions, the optimal choice 
is (Lilly 1992): 


( MijLij ) ~ f 2 ^l2~~ ' 

c s = ~ 4 , Mij = 2 A \S\Sij - A \S\Sij , 

\MijMijj L 


( 2 . 11 ) 
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where (• • • ) denotes spatial averaging. Thus, although model (2.7) does not explicitly 
incorporate the magnetic field, the dynamic constant C s should adjust to an appropriate 
value if the scale-similarity hypothesis is at all justified. 


3. Numerical experiments 

To assess the dynamic Smagorinsky model in the present context, we have built a 
set of DNS databases. Since we restrict our attention to a cubic domain with periodic 
boundary conditions, a spectral (dealiased) code has been used. 

The velocity field can then be initialized in Fourier space. The initial mode amplitudes 
are computed to match the spectra of the Comte-Bellot & Corrsin (1971) experiment at 
stage 1 (see Rogallo (1981) for more information on the procedure). Phases are initially 
random and the flow is left to freely decay until the skewness of the velocity derivative 
reaches a quasi constant value of S = —0.4. At this time (hereafter referred to as to) 
the flow is considered “physical” and used as the initial condition for all of our runs. All 
the DNS simulations are done in a (27t) 3 box using a resolution of 512 3 Fourier modes 
and the viscosity is set to v = 0.006. Other relevant quantities measured at t = to are 
summarized here: 



380 (integral Reynolds number), 



3 u 2 

Teddy = ^ ~ = 


= 84.1 (microscale Reynolds number), 
0.238 (eddy turnover time). 


(3.1) 

(3.2) 

(3.3) 


Three test cases have been considered. The first one corresponds to a decaying flow 
without the addition of any applied external magnetic field. The other two numerical 
experiments are distinguished by their characteristic interaction numbers at to- The in- 
teraction number N (also known as the Stuart number) is defined as follows: 


N = 


(Bt xt ) 2 L 
rj u 


(3.4) 


N measures the relative strengths of the magnetic damping term and the non-linear term 
in (2.2). Two cases are examined here, N = 1 and N = 10. 

To illustrate the effect of the external magnetic field, the time history of the kinetic 
energy density E = y f dx^Ui(x)ui(x) of the flow for the three cases is presented in Fig. 
1. As in all the following figures, time has been normalized by the initial eddy turnover 
time and the non-dimensional time is denoted by t*. 


4. LES results 

In order to assess the LES method, some DNS snapshots of the flow field have been 
truncated from the 512 3 resolution to a 32 3 resolution (using a sharp Fourier cut-off). 
The initial condition for the LES runs has been obtained in the same way by truncating 
the DNS field at t = t 0 . As can be easily seen from (2.11), the only free parameter in the 

subgrid-scale model is the ratio A/ A. As is standard practice in LES of non-conductive 
flows, the value of 2 is adopted here. 
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Figure 1. Time history of the kinetic energy density E. The solid line represents the flow 
decay without applied magnetic field, the dashed curve corresponds to the case N = 1 and the 
dasli-dot curve corresponds to the case N — 10. The dotted line represents the time to at which 
the magnetic field is switched on. 





t* 


t* 


t* 


Figure 2. Time history of the resolved kinetic energy density: LES vs. filtered DNS. The solid 
lines represent the predictions of the LES. The diamonds represent the filtered kinetic energy 
density obtained by truncating the DNS fields to a 32 3 resolution. The dotted lines represent 
the predictions of the unresolved DNS (performed on the 32 3 mesh) without any subgrid-scale 
model. 


4.1. Kinetic Energy 

Figure 2 represents the time evolution of the resolved kinetic energy density Er = 
vf dx gUj(x)uj(x) predicted by the LES and compared to the filtered DNS. On each 
plot, a third curve representing a simulation (referred to as ‘unresolved DNS’) on the 
32 3 mesh without SGS modeling is added to stress the action of the subgrid-scale model. 
The case B ext = 0 serves as a benchmark to check that in the case of non-conductive 
flows, our LES code behaves as expected. In both the cases N = 1 and N = 10, the LES 
performs remarkably well. In the case N = 1 the difference between LES and unresolved 
DNS is very clear. In the case N = 10 and for this diagnostic, the unresolved DNS does 
not depart significantly from the filtered DNS and LES. 

4.2. Energy spectra 

The energy spectrum constitutes a finer diagnostic than the kinetic energy density since 
it retains information about the repartition of energy among the different scales of the 
flow. Figures 3 and 4 contain respectively kinetic energy spectra for the cases TV = 1 
and TV = 10 at several instants in the simulation. Agreement between LES and DNS is 
again excellent. The unresolved DNS exhibits the usual pile up of energy at high wave 
numbers that results from lack of resolution. Thus, although in terms of global energy 
the unresolved DNS does perform quite well in the TV = 10 case, the spectral properties 
of the flow are poorly predicted and the LES does a much better job. 
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IGURE 3. Three-dimensional kinetic energy spectra for the case N — 1. The solid line represents 
re filtered DNS, the dashed line represents the LES and the dotted line corresponds to the 
nresolved DNS. The times at which the spectra are calculated are indicated in the plots. 
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Figure 4. Three-dimensional kinetic energy spectra for the case N = 10. The solid line repre- 
sents the filtered DNS, the dashed line represents the LES and the dotted line corresponds to 
the unresolved DNS. The times at which the spectra are calculated are indicated in the plots. 


4.3. Dissipation rate budget 

The evolution of the resolved kinetic energy density E B is governed by, 


dEn 

dt 


1 

V 


dx.(e v + e B + e sgs ) 


where, 


— 2 u Sjj S jj . 


= 


{B. 


ext\2 


iA d z d z 


c sgs 


— TijSij. 


(4.1) 


(4.2) 


From Fig. 2 it is clear that the LES predicts the total dissipation rate very well. It is 
however instructive to know how the different contributions (4.2) are reproduced sepa- 
rately. 

In Fig. 5 we present the time evolution of the resolved viscous dissipation rate e„ for the 
LES, the filtered DNS and the unresolved DNS. The figures show that the LES reproduces 
this diagnostic very well, whereas the unresolved DNS systematically overestimates it. 
This results from the pile-up of energy near the LES cut-off where most of the resolved 
viscous dissipation occurs. 

Figure 6 represents the time evolution of the resolved magnetic dissipation e B . Again 
the LES predictions match the filtered DNS results very well. The differences with the 
unresolved DNS are less severe than for the viscous dissipation. This is to be expected 
since the magnetic dissipation occurs at every scale in the flow and its overall intensity 
is thus less contaminated by the pile-up of energy occurring near the LES cut-off for the 
unresolved DNS. It is also interesting to note that in the case N = 10, the magnetic 
dissipation falls very rapidly with time at the early stages of the decay. This happens 
because in that case the flow quickly becomes fairly independent of the 2 -direction which 
is parallel to the magnetic field (this is illustrated further in Section 4.4). 

Finally, in Fig. 7 the evolution of the subgrid-scale transfer rate e sgs is presented. For 
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Figure 5. Resolved viscous dissipation rate e„ for the cases TV = 1 (left) and TV = 10 (right). 
In each plot, the solid line represents the LES, the diamonds represent the filtered DNS and the 
dotted line represents the unresolved DNS. 




Figure 6. Resolved magnetic dissipation rate es for the cases TV = 1 (left) and TV = 10 (right). 
In each plot, the solid line represents the LES, the diamonds represent the filtered DNS and the 
dotted line represents the unresolved DNS. 


the LES, e sgs is obtained by substituting the dynamic Smagorsinky in place of r,j in 
(4.2). DNS results are obtained by computing the exact fy = UiUj — uiUj from the DNS 
data fields. Obviously, no (non-zero) results are available for the unresolved DNS. From 
the plots we see that the subgrid-scale transfer rate is initially overestimated in the LES. 
This is not a surprise since the dynamic procedure usually needs a little time to settle. 
Soon after this short transient time, the agreement between LES and DNS predictions 
is very good. 


4.4. Flow structures 

It is well known that the extra damping term present in the quasi-static approximation 
(2.2) leads to a progressive suppression of spatial variations in the flow along the direction 
of the magnetic field. It is thus important to assess how well the LES is able to reproduce 
this feature. To that end, we have plotted in Figures 8 and 9 the contours of the kinetic 
energy density at three different times respectively for the filtered DNS and the LES 
(only the case TV = 10 is shown because the effect is more pronounced for strong applied 
magnetic fields). Since the turbulence is decaying with time, the colormap had to be 
rescaled for each of the single plots. However, the same colormap for the filtered DNS 
and LES has been used in the plots corresponding to the same times. As is obvious from 
the plots, the LES reproduces the filtered DNS structures very well. The unresolved DNS 
also captures the large-scales structures of the flow (no plots shown) but due to the pile 
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Figure 7. Subgrid-scale transfer rate e sgs for the cases N = 1 (left) and N = 10 (right). In each 
plot, the solid line represents results obtained from the LES while the diamonds are obtained 
from the DNS data. 


t* = 0.285 t* = 1.02 t* = 2.68 

Figure 8. Contours of the kinetic energy obtained from the filtered DNS ( N = 10 case). The 
different times at which the contours are calculated are indicated under the plots and are the 
same ones as those used in Fig. 9. 


up of energy it predicts a significant amount of small scales; the corresponding contour 
plots contain some excessive ‘noise’ on top of the large-scales structures. 


4.5. Dynamic Constant 

As recalled in Section 2, the dynamic procedure is designed to optimize the scaling 
constant C s present in the Smagorinsky subgrid-scale model. We also mentioned that 
when an external magnetic field is present, the value of C s should decrease since non- 
linear transfers are reduced and the effect of the subgrid-scale model should be damped 


t* = 0.285 t* = 1.02 t* = 2.68 

Figure 9. Contours of the kinetic energy obtained from the LES ( N = 10 case). The different 
times at which the contours are calculated are indicated under the plots and are the same ones 
as those used in Fig. 8. 
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Figure 10. Time history of the dynamic Smagorinsky constant. The solid line represents the 
case with no magnetic field, the dashed line corresponds to the case N = 1 and the case N =10 
is represented by the dotted line. 

2 

accordingly. In Fig. 10 the evolution of C s A is plotted. From the figure, it is clear that 
the expected behavior is observed. As the interaction number is increased, the dynamic 
procedure automatically adapts the value of the Smagorinsky constant to decrease the 
effect of the subgrid-scale model. 

5. Conclusions 

In this article we have shown that the dynamic Smagorinsky model can be used to 
perform large-eddy simulations of flows subject to an applied external magnetic field at 
low magnetic Reynolds number, R m . Although this subgrid-scale model was not designed 
for this application, its behavior is excellent owing to its adaptation to the flow and the 
applied magnetic field through the dynamic procedure. The model can be considered 
robust since it works equally well for interaction (or Stuart) numbers ranging from N = 0 
(no magnetic field) to N = 10 (for which the flow becomes nearly two dimensional). In 
the future, the same model will be tested in more complex geometries. 
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